#include "generalization.h"
#include <iostream>
#include <device_launch_parameters.h>
#include <thrust/device_vector.h>
#include <thrust/sequence.h>
#include <thrust/sort.h>
#include <fstream>
#include <thrust/host_vector.h>
#include <thrust/scan.h>
#include <fstream>
#include <stdlib.h>     
#include <time.h> 
#include <math.h>
#include <set>
#include <sstream>
#include <device_launch_parameters.h>
#include <algorithm>
#include <cuda_profiler_api.h>
using namespace std;
using namespace thrust;


vector<double> exec_time;

// compare upper bound of Result object
struct Ubcomapre
{
	__host__ __device__ bool operator()(const Result &lhs, const Result &rhs) const{
		if(lhs.query < rhs.query) return true;
		else if(lhs.query > rhs.query) return false;
		if(lhs.ub == rhs.ub) return lhs.lb < rhs.lb;
		else return lhs.ub < rhs.ub;
	}
};

// compare lower bound of Result object
struct Lbcompare
{
	__host__ __device__ bool operator()(const Result &lhs, const Result &rhs) const{
		if(lhs.query < rhs.query) return true;
		else if(lhs.query > rhs.query) return false;
		return lhs.lb < rhs.lb;
	}
};

bool compare_result(const Result& r1, const Result& r2){
	return r1.query < r2.query || (r1.query == r2.query && r1.feature_id < r2.feature_id);
}

// distance function for values on the same dimension
template <class T>
__device__ T distance_func(T& a, T& b){
	switch(D_MAX_FEATURE_ID){
	case 1:
	{
		T diff = a - b;
		return diff*diff;
	}
	default:
	{
		T diff = a - b;
		return diff*diff;
	}
	}
}

// allocate GPU memory
template <class T>
void init_gpu_memory(host_vector<T>& hvec, device_vector<T>** dvec){
	*dvec = new device_vector<T>(hvec.size());
	**dvec = hvec;
}

// get a random number
double get_uniform_rand(){
	return (double)(rand() % 100000) / 100000;
}

// computing the counting while saving the position
__global__ void compute_mapping_saving_pos(
	int* query_keyword,
	InvlistEnt* invert_list, int* invert_list_idx,
	QueryFeatureEnt* query_feature,
	int* query_search_dims,int* query_search_dims_idx,
	int N, // number of elements processd this round
	QueryPos* query_pos, // remeber the counting position
	int lb, // lower bound of the search, relative bounds 
	int ub // upper bound of the search, relative bounds
	){

	int bid = blockIdx.x; // correspond to query
	int tid = threadIdx.x;
	int block_size = blockDim.x;

	int search_dim_start = (bid == 0) ? 0 : query_search_dims_idx[bid-1];
	int search_dim_end = query_search_dims_idx[bid];

	int start_idx = bid * MAX_DIM;

	// loading the position
	__shared__ int2 pos[MAX_DIM];
	if(tid < MAX_DIM){
		pos[tid] = query_pos[bid].pos[tid];
	}

	__syncthreads();

	for(int iter = search_dim_start; iter < search_dim_end; iter++){
		int dim = query_search_dims[iter];
		int down_pos = pos[dim].x;
		int up_pos = pos[dim].y;

		int keyword = query_keyword[start_idx+dim];
		int y_dim_value = keyword % MAX_DIM_VALUE;

		if(down_pos == 0 && up_pos == 0){
			int invert_list_start = keyword == 0 ? 0 : invert_list_idx[keyword-1];
			int invert_list_end = invert_list_idx[keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					int target_idx = bid*D_MAX_FEATURE_ID+inv_ent;
					query_feature[target_idx].count += 1;
				}
			}
		}


		// going downward
		int down_compute = N;
		int down_value = y_dim_value - down_pos; // move the position that last iteartion possessed.

		while(down_value - 1 >= 0 && down_compute > 0 && down_value > y_dim_value - lb){
			down_value--;
			int down_keyword = down_value + dim * MAX_DIM_VALUE;

			int invert_list_start = down_keyword == 0 ? 0 : invert_list_idx[down_keyword-1];
			int invert_list_end = invert_list_idx[down_keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			// calcuate the distance from the current point to the query along this dimension
			int dist_sqr = distance_func(y_dim_value,down_value);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					int target_idx = bid*D_MAX_FEATURE_ID+inv_ent;
					query_feature[target_idx].count += 1;
					query_feature[target_idx].ACD += dist_sqr;
				}
			}

			down_compute -= invert_list_size;
		}

		// going upward
		int up_compute = N;
		int up_value = y_dim_value + up_pos; // move the position that last iteartion possessed.

		while(up_value + 1 < MAX_DIM_VALUE && up_compute > 0 && up_value < y_dim_value + ub){
			up_value++;
			int up_keyword = up_value + dim * MAX_DIM_VALUE;

			int invert_list_start = up_keyword == 0 ? 0 : invert_list_idx[up_keyword-1];
			int invert_list_end = invert_list_idx[up_keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			// calcuate the distance from the current point to the query along this dimension
			int dist_sqr = distance_func(up_value,y_dim_value);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					int target_idx = bid*D_MAX_FEATURE_ID+inv_ent;
					query_feature[target_idx].count += 1;
					query_feature[target_idx].ACD += dist_sqr;
				}
			}

			up_compute -= invert_list_size;
		}

		__syncthreads();
		if(tid == 0){
			pos[dim].x =  y_dim_value - down_value;
			pos[dim].y = up_value - y_dim_value;
		}
	}

	__syncthreads();
	if(tid < MAX_DIM){
		query_pos[bid].pos[tid] = pos[tid];
	}
}

__global__ void prefix_count(
	QueryFeatureEnt* query_feature,
	int* count_vec, int* blk_count_vec, float threshold,int round_num){

	int bid = blockIdx.x;
	int tid = threadIdx.x;
	int oid = bid * blockDim.x + tid;

	int block_start = bid * D_MAX_FEATURE_ID;
	int block_end = block_start + D_MAX_FEATURE_ID;
	int round = D_MAX_FEATURE_ID / blockDim.x + (D_MAX_FEATURE_ID % blockDim.x != 0);


	__shared__ int shared_count;
	if(tid == 0) shared_count = 0;
	__syncthreads();

	int count = 0;
	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + tid;
		if(idx < block_end){
			if(query_feature[idx].count > threshold){
				count++;
				atomicAdd(&shared_count,1);
			}
		}
	}
	count_vec[oid] = count;
	__syncthreads();

	if(tid == 0){
		blk_count_vec[bid] = shared_count;
	}
}

__device__ void blk_sort_shared_memory(int* shared_array, int* temp_array, int size){
	int rounds = ceil(log2((double)size));

	int base = 1;
	for(int i = 1; i <= rounds; i++){
		
		if(threadIdx.x % (2*base) == 0){
			int start_idx_x = threadIdx.x;
			int end_idx_x = start_idx_x + base;
			end_idx_x = end_idx_x < size ? end_idx_x : size;

			int start_idx_y = end_idx_x;
			int end_idx_y = start_idx_y + base;
			end_idx_y = end_idx_y < size ? end_idx_y : size;

			int x_pointer = start_idx_x;
			int y_pointer = start_idx_y;

			int output_ptr = x_pointer;
			while(x_pointer < end_idx_x || y_pointer < end_idx_y){
				if(x_pointer >= end_idx_x)
					temp_array[output_ptr++] = shared_array[y_pointer++];
				else if(y_pointer >= end_idx_y)
					temp_array[output_ptr++] = shared_array[x_pointer++];
				else if(shared_array[x_pointer] < shared_array[y_pointer])
					temp_array[output_ptr++] = shared_array[y_pointer++];
				else
					temp_array[output_ptr++] = shared_array[x_pointer++];
			}			
		}
		__syncthreads();
		if(threadIdx.x < size){
			shared_array[threadIdx.x] = temp_array[threadIdx.x];
		}
		__syncthreads();
		base *= 2;
	}
}

__global__ void output_result(
	QueryFeatureEnt* query_feature,
	QueryPos* query_pos, // remeber the counting position
	int* query_keyword, // the information of query and it's value for each dimension
	int* ending_idx,float threshold,
	Result* result_vec, int round_num){

	int pos = blockIdx.x * blockDim.x + threadIdx.x;
	
	int start_idx = pos == 0 ? 0 : ending_idx[pos-1];

	int block_start = blockIdx.x * D_MAX_FEATURE_ID;
	int block_end = block_start + D_MAX_FEATURE_ID;
	int round = D_MAX_FEATURE_ID / blockDim.x + (D_MAX_FEATURE_ID % blockDim.x != 0);

	__shared__ int queryLowerBound[MAX_DIM];
	__shared__ int queryUpperBound[MAX_DIM];
	__shared__ int temp_shared[MAX_DIM];

	if(threadIdx.x < MAX_DIM){
		int2 q_pos = query_pos[blockIdx.x].pos[threadIdx.x];
		int max_lb = max(q_pos.x,q_pos.y);
		queryLowerBound[threadIdx.x] = max_lb * max_lb;

		int query_start_idx = blockIdx.x * MAX_DIM;
		int query_dim_value = (query_keyword[query_start_idx + threadIdx.x]) % MAX_DIM_VALUE;
		int max_ub =  max(query_dim_value,MAX_DIM_VALUE-query_dim_value);
		queryUpperBound[threadIdx.x] = max_ub * max_ub;
	}

	__syncthreads();

	// sort the lower and up bound
	blk_sort_shared_memory(queryLowerBound,temp_shared,MAX_DIM);
	blk_sort_shared_memory(queryUpperBound,temp_shared,MAX_DIM);
	
	// calculate lower bound for this query, total distance for current retrieve point
	int lower_bound_sum = 0;
	for(int i = 0; i < MAX_DIM; i++) lower_bound_sum += queryLowerBound[i];

	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + threadIdx.x;
		if(idx < block_end){
			float count = query_feature[idx].count;
			if(count > threshold){
				Result new_result;
				new_result.query = blockIdx.x;
				new_result.feature_id = i * blockDim.x + threadIdx.x;
				new_result.count = count;
				
				int ACD = query_feature[idx].ACD;

				// DL(f) = DL(us) - sum(max count number of dt) + ACD(f)
				int lb = lower_bound_sum;
				for(int m = 0; m < (int)count; m++){
					lb -= queryLowerBound[m];
				}
				lb += ACD;

				// DU(f) = ACD(f) + sum(max 128 - c of upper bound)
				int ub = ACD;
				for(int m = 0; m < MAX_DIM - (int)count; m++){
					ub += queryUpperBound[m];
				}

				new_result.lb = lb;
				new_result.ub = ub;

				result_vec[start_idx] = new_result;
				start_idx++;
			}
		}
	}
}

__global__ void terminate_check(Result* ub_sorted,Result* lb_sorted, int* end_idx, int* output, int K, int round_num){
	int tid = threadIdx.x;
	int bid = blockIdx.x;

	int blk_start_idx = bid == 0 ? 0 : end_idx[bid-1];
	int blk_end_idx = end_idx[bid];

	/*if(tid == 0 && bid == 26){
		printf("reached here: bid: %d with block start idx %d and end idx %d \n",bid,blk_start_idx,blk_end_idx);
	}*/

	if(blk_end_idx - blk_start_idx <= K){
		return;
	}

	blk_end_idx = blk_start_idx + K + 1;

	int round = (K+1) / blockDim.x + ((K+1) % blockDim.x != 0);

	__shared__ int termination_check;
	if(tid == 0) termination_check = 0;
	__syncthreads();

	for(int i = 0; i < round; i++){
		int idx = blk_start_idx + i * blockDim.x + tid;
		if(idx < blk_end_idx){
			Result lb_elem = lb_sorted[idx];
			bool duplicated = false; // if element in the lb sorted array existed in the ub sorted array
			for(int j = 0; j < K; j++){
				Result ub_elem = ub_sorted[blk_start_idx + j];
				if(lb_elem.feature_id == ub_elem.feature_id){
					duplicated = true;
					break;
				}
			}
			Result ub_elem = ub_sorted[blk_start_idx + K - 1];
			bool lb_greater_ub = (lb_elem.lb >= ub_elem.ub);

			if(!duplicated && lb_greater_ub){
				//printf("reached here tid %d and bid %d\n",tid,bid);
				atomicAdd(&termination_check,1);
			}
		}
	}
	__syncthreads();
	if(tid == 0){
		if(termination_check != 0){
			output[bid] = 0;
		}
	}
}

__global__ void extract_topK(Result* ub_sorted,int* end_idx,Result* output, int K){
	int tid = threadIdx.x;
	int bid = blockIdx.x;
	int oid = bid * K + tid;

	int blk_start_idx = bid == 0 ? 0 : end_idx[bid-1];
	int blk_end_idx = blk_start_idx + K;

	int round = K / blockDim.x + (K % blockDim.x != 0);

	for(int i = 0; i < round; i++){
		int idx = blk_start_idx + i * blockDim.x + tid;
		if(idx < blk_end_idx){
			output[oid] = ub_sorted[idx];
			oid += blockDim.x;
		}
	}
}

QueryManager::QueryManager(string inv_list_dir,int distance_func){
	if(inv_list_dir.empty()) rand_inv_list();
	else get_invert_list_from_file(inv_list_dir);

	cudaMemcpyToSymbol(D_MAX_FEATURE_ID,&this->MAX_FEATURE_ID,sizeof(int));
	cudaMemcpyToSymbol(D_DIST_FUNCTION_TYPE,&distance_func,sizeof(int));

	query_initialized = false;
	d_query_keyword = NULL;
	d_query_feature = NULL;
	d_query_pos = NULL;
	d_valid_query = NULL;
	d_query_search_dims = NULL;
	d_query_search_dims_idx = NULL;
}

void QueryManager::init_new_query(string query_file_path,int num_of_query,int N,int lb,int ub,float threshold, int K){
	clear_query_memory();
	
	cout << "start initilizing query" << endl;

	// read in read queries
	if(query_file_path.empty())rand_query(num_of_query);
	else{
		cout << "this feature is not supported now since query search dimension is not implemented" << endl;
		exit(1);
		get_query_from_file(query_file_path);
	}

	// checking memory size
	if((double)MAX_FEATURE_ID*d_query_pos->size()*sizeof(QueryFeatureEnt) > 5.0 * pow(2.0,30.0)){
		cout << "too much memory is used for query feature mapping" << endl;
		exit(1);
	}

	/*double total_bytes = (double)d_invert_list->size() * sizeof(InvlistEnt) + 
	d_invert_list_idx->size() * sizeof(int) + 
	d_query_keyword->size() * sizeof(int) + 
	d_query_pos->size() * sizeof(QueryPos) + 
	MAX_FEATURE_ID*d_query_pos->size()*sizeof(int);*/

	d_query_feature = new device_vector<QueryFeatureEnt>((unsigned int)MAX_FEATURE_ID*d_query_pos->size());
	d_valid_query = new device_vector<int>(d_query_pos->size(),1);

	num_of_rounds = 0;
	this->N = N;
	this->lb = lb;
	this->ub = ub;
	this->threshold = threshold;
	this->K = K;

	this->query_initialized = true;
	cout << "finish initialization" << endl;
}

void QueryManager::clear_query_memory(){
	if(d_query_keyword != NULL) delete d_query_keyword;
	if(d_query_feature != NULL) delete d_query_feature;
	if(d_query_pos != NULL) delete d_query_pos;
	if(d_valid_query != NULL) delete d_valid_query;
	if(d_query_search_dims != NULL) delete d_query_search_dims;
	if(d_query_search_dims_idx != NULL) delete d_query_search_dims_idx;
}

QueryManager::~QueryManager(){
	delete d_invert_list;
	delete d_invert_list_idx;
	clear_query_memory();
}

void QueryManager::rand_inv_list(){	
	MAX_FEATURE_ID = 1000000;
	vector<vector<int> > keyword_feature(MAX_DIM*MAX_DIM_VALUE); 
	for(int i = 0; i < MAX_FEATURE_ID; i++){
		for(int x = 0; x < MAX_DIM;  x++){
			int keyword = x * MAX_DIM_VALUE + rand() % MAX_DIM_VALUE;
			keyword_feature[keyword].push_back(i);
		}
	}

	host_vector<InvlistEnt> h_invert_list;
	host_vector<int> h_invert_list_idx;
	int count = 0;
	// construct inverted list
	for(int i = 0; i < MAX_DIM*MAX_DIM_VALUE; i++){
		for(int j = 0; j < keyword_feature[i].size(); j++){
			InvlistEnt new_ent;
			new_ent = keyword_feature[i][j];
			h_invert_list.push_back(new_ent);
		}
		count += keyword_feature[i].size();
		h_invert_list_idx.push_back(count);
	}
	
	// transfer to gpu
	init_gpu_memory(h_invert_list,&d_invert_list);
	init_gpu_memory(h_invert_list_idx,&d_invert_list_idx);
}

void QueryManager::get_query_from_file(string query_file_name){
	ifstream query_file(query_file_name.c_str());
	if(query_file.is_open()){
		cout << "query file cannot found" << endl;
		exit(1);
	}

	host_vector<int> h_query_keyword;
	int num_query = 0;
	while(!query_file.eof()){
		string input_line;
		getline(query_file,input_line);
		stringstream ss_test(input_line);

		string test;
		ss_test >> test;
		if(test.empty()) continue;

		num_query++;
		stringstream ss(input_line);

		for(int i = 0; i < MAX_DIM; i++){
			int keyword = -1;
			ss >> keyword;
			h_query_keyword.push_back(keyword);
		}
	}

	d_query_pos = new device_vector<QueryPos>(num_query);
	init_gpu_memory(h_query_keyword,&d_query_keyword);
	cout << "finish reading query from file" << endl;
}	

void QueryManager::rand_query(int query_num){

	host_vector<int> h_query_keyword;
	// construct feature to keyword
	for(int i = 0; i < query_num; i++){
		for(int x = 0; x < MAX_DIM; x++){
			int new_keyword = x * MAX_DIM_VALUE + rand() % MAX_DIM_VALUE;
			h_query_keyword.push_back(new_keyword);
		}
	}

	d_query_pos = new device_vector<QueryPos>(query_num);

	init_gpu_memory(h_query_keyword,&d_query_keyword);

	// construct search dimension for query
	vector<set<int> > search_dims(query_num);
	host_vector<int> h_query_search_dims;
	host_vector<int> h_query_search_dims_idx;
	int count = 0;
	for(int i = 0; i < query_num; i++){
		/*int num_dims = rand() % MAX_DIM;
		while(search_dims[i].size() != num_dims){
			search_dims[i].insert(rand() % MAX_DIM);
		}
		count += num_dims;
		h_query_search_dims_idx.push_back(count);
		for(set<int>::iterator it = search_dims[i].begin(); it != search_dims[i].end(); it++)
			h_query_search_dims.push_back(*it);*/
		count += MAX_DIM;
		h_query_search_dims_idx.push_back(count);
		for(int j = 0; j < MAX_DIM; j++) h_query_search_dims.push_back(j);
	}

	init_gpu_memory(h_query_search_dims,&d_query_search_dims);
	init_gpu_memory(h_query_search_dims_idx,&d_query_search_dims_idx);
}

void QueryManager::get_invert_list_from_file(string data_root_dir){

	string sorted_file_name(data_root_dir);
	sorted_file_name.append("sorted.txt");
	ifstream sorted_file(sorted_file_name.c_str());
	if(!sorted_file.is_open()){
		cout << "sorted index file not found" << endl;
		exit(1);
	}

	cout << "start to read invert list" << endl;

	host_vector<InvlistEnt> h_inv_list;
	host_vector<int> h_inv_list_idx;
	int count = 0;
	int keyword_idx = 0;

	MAX_FEATURE_ID = -1;

	while(!sorted_file.eof()){
		string input_line;
		getline(sorted_file,input_line);
		stringstream ss(input_line);

		string name;
		int file_idx;
		int file_len;

		ss >> name;
		if(name.empty()) continue;
		ss >> file_len;

		file_idx = atoi(name.c_str());
		string file_name(data_root_dir);
		file_name.append("sift_small/");
		file_name.append(name.c_str());

		for(int i = keyword_idx; i < file_idx; i++){
			h_inv_list_idx.push_back(count);
		}

		keyword_idx = file_idx;
		count += file_len;

		ifstream keyword_file(file_name.c_str());
		if(!keyword_file.is_open()){
			cout << "invert list file not found" << endl;
			exit(1);
		}

		for(int i = 0; i < file_len; i++){
			InvlistEnt new_ent;
			keyword_file >> new_ent;
			h_inv_list.push_back(new_ent);
			if(new_ent > MAX_FEATURE_ID){
				MAX_FEATURE_ID = new_ent;
			}
		}

		keyword_file.close();

	} 
	sorted_file.close();
	cout << "finish reading invert list" << endl;

	MAX_FEATURE_ID++;

	for(int i = keyword_idx; i < MAX_DIM * MAX_DIM_VALUE; i++){
		h_inv_list_idx.push_back(count);
	}

	init_gpu_memory(h_inv_list,&d_invert_list);
	init_gpu_memory(h_inv_list_idx,&d_invert_list_idx);
}


bool QueryManager::run_query(vector<Result>& result){
	if(!this->query_initialized){
		cout << "query has not been initialized" << endl;
		exit(1);
	}
	if(THREAD_PER_BLK < MAX_DIM){
		cout << "thread per block must be larger than max dimension of x" << endl;
		exit(1);
	}
	if(d_query_pos->size() > 65535){
		cout << "number of query should not exceed 65535: need to modify the code" << endl;
		exit(1);
	}
	if(THREAD_PER_BLK > 1024){
		cout << "THREAD_PER_BLK should not exceed 1024: need to modify the code" << endl;
		exit(1);
	}

	cudaProfilerStart();
	//cout << "started counting" << endl;

	cudaEvent_t start, stop;
	float elapsedTime;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);

	cudaEventRecord(start, 0);

	compute_mapping_saving_pos<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_query_keyword->data()),
		raw_pointer_cast(d_invert_list->data()),
		raw_pointer_cast(d_invert_list_idx->data()),
		raw_pointer_cast(d_query_feature->data()),
		raw_pointer_cast(d_query_search_dims->data()),
		raw_pointer_cast(d_query_search_dims_idx->data()),
		N,raw_pointer_cast(d_query_pos->data()),lb,ub);

	device_vector<int> threshold_count(d_query_pos->size()*THREAD_PER_BLK,0);
	device_vector<int> query_result_count(d_query_pos->size(),0);
	//cout << "started prefix counting" << endl;

	prefix_count<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_query_feature->data()),
		raw_pointer_cast(threshold_count.data()),
		raw_pointer_cast(query_result_count.data()),
		threshold,num_of_rounds);

	//if(num_of_rounds == 18){
		/*thrust::host_vector<QueryFeatureEnt> dvec(d_query_feature->begin(),d_query_feature->end());
		cout << "count at query 26 and feature 0:" << dvec[26*MAX_FEATURE_ID].count << endl;
		cout << "%%%%%%%%%%%%%%%%" << endl;*/
	//}
		/*thrust::host_vector<QueryPos> ddddvec(d_query_pos->begin(),d_query_pos->end());
		for(int i = 0; i < 128; i++){
			cout << "(" << ddddvec[26].pos[i].x << "," << ddddvec[26].pos[i].y << ") ";
		}
		cout << endl;
		*/
	//	thrust::host_vector<InvlistEnt> dinv(d_invert_list->begin(),d_invert_list->end());
		/*thrust::host_vector<int> dinvidx(d_invert_list_idx->begin(),d_invert_list_idx->end());
		cout << "invert list size: " << dinvidx.size() << endl;*/

		/*for(int i = 0; i < 128; i++){
			int count = 0;
			for(int j = 0; j < 256; j++){
				int key = i*256+j;
				int start_idx = (key == 0) ? 0 : dinvidx[key-1];
				int end_idx = dinvidx[key];
				for(int k = start_idx; k < end_idx; k++){
					if(dinv[k].feature_id == 0) count++;
				}
			}
			if(count != 1) cout << " dim: " << i << " has zero: " << count;
		}*/
		
	thrust::inclusive_scan(threshold_count.begin(),threshold_count.end(),threshold_count.begin()); // per thread inclusive scan
	thrust::inclusive_scan(query_result_count.begin(),query_result_count.end(),query_result_count.begin()); // per block inclusive scan

	device_vector<Result> d_result_ub_sorted(threshold_count[threshold_count.size()-1]);
	
	//cout << "start output result" << endl;

	output_result<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_query_feature->data()),
		raw_pointer_cast(d_query_pos->data()),
		raw_pointer_cast(d_query_keyword->data()),
		raw_pointer_cast(threshold_count.data()),
		threshold,
		raw_pointer_cast(d_result_ub_sorted.data()),num_of_rounds);
	
	device_vector<Result> d_result_lb_sorted(d_result_ub_sorted.begin(),d_result_ub_sorted.end());

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[0] += elapsedTime;
	
	cudaEventRecord(start, 0);

	thrust::sort(d_result_ub_sorted.begin(),d_result_ub_sorted.end(),Ubcomapre());
	thrust::sort(d_result_lb_sorted.begin(),d_result_lb_sorted.end(),Lbcompare());

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[1] += elapsedTime;

	/*double time = (double)(clock() - start) / CLOCKS_PER_SEC;
	cout << "sort timing: " << time << endl;*/

	cudaEventRecord(start, 0);
	terminate_check<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_result_ub_sorted.data()),
		raw_pointer_cast(d_result_lb_sorted.data()),
		raw_pointer_cast(query_result_count.data()),
		raw_pointer_cast(d_valid_query->data()),K,num_of_rounds);  

	int terminate_sum = thrust::reduce(d_valid_query->begin(),d_valid_query->end());

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[2] += elapsedTime;
	cudaEventDestroy(start);
	cudaEventDestroy(stop);

	cout << "number of total result: " << d_result_ub_sorted.size() << " with " << terminate_sum << " queries still running in rounds " << num_of_rounds << endl;

	cudaProfilerStop();
	num_of_rounds++;
	
	device_vector<Result> d_result(d_query_pos->size()*K);
	extract_topK<<<d_query_pos->size(),THREAD_PER_BLK>>>(
		raw_pointer_cast(d_result_ub_sorted.data()),
		raw_pointer_cast(query_result_count.data()),
		raw_pointer_cast(d_result.data()),K);

	//return result
	result.resize(d_result.size());
	cudaMemcpy(&result[0],raw_pointer_cast(d_result.data()),d_result.size()*sizeof(Result),cudaMemcpyDeviceToHost);
	
	return terminate_sum == 0;
}


void QueryManager::print_problem(){
	cout << "query to keyword relationship" << endl;
	for(int i = 0; i < d_query_pos->size(); i++){
		cout << "Q " << i << ":";
		for(int j = 0; j < MAX_DIM; j++){
			int keyword = (*d_query_keyword)[i*MAX_DIM+j];
			int x = keyword / MAX_DIM_VALUE;
			int y = keyword % MAX_DIM_VALUE;
			cout << "(" << x << "," << y << ");";
		}
		cout << endl;
	}
	cout << "inverted list" << endl;
	for(int x = 0; x < MAX_DIM; x++){
		for(int y = 0; y < MAX_DIM_VALUE; y++){
			

			int keyword = x*MAX_DIM_VALUE+y;
			int start_idx = (keyword == 0) ? 0 : (*d_invert_list_idx)[keyword-1];
			int end_idx = (*d_invert_list_idx)[keyword];

			if(start_idx == end_idx) continue;

			cout << "(" << x << "," << y << "):";
			for(int i = start_idx; i < end_idx; i++){
				InvlistEnt ent = (*d_invert_list)[i];
				cout << ent << " ";
			}
			cout << endl;
		}
	}
}

void QueryManager::print_final_result(vector<Result>& result){
	cout << "result size is: " << result.size() << endl;
	std::sort(result.begin(),result.end(),compare_result);
	int qindex = -1;
	for(int i = 0; i < result.size(); i++){
		Result new_result = result[i];
		if(new_result.query != qindex){
			cout << endl;
			cout << new_result.query << ":" << "(" << new_result.feature_id << "," << new_result.count << ") ";
			qindex = new_result.query;
		}else{
			cout << "(" << new_result.feature_id << "," << new_result.count << ") ";
		}
	}
}

int main(){
	
	srand(1);

	vector<Result> result;
	int N = 1000;
	int num_of_query = 1000;
	int lb = 256;
	int ub = 256;
	float threshold = 3;
	int K = 3;

	QueryManager query_handle("../data/sift_exp_data/");
	
	query_handle.init_new_query("",num_of_query,N,lb,ub,threshold,K);
	cout << "finish init" << endl;

	exec_time.resize(3,0);

	double total_time = 0;
	bool finished = false;
	while(!finished){
		
		clock_t start = clock();
		finished = query_handle.run_query(result);
		cudaDeviceSynchronize();
		double time = (double)(clock() - start) / CLOCKS_PER_SEC;
		total_time += time;
		cout << "iteration takes time: " << time << endl;
	}

	query_handle.init_new_query("",num_of_query,N,lb,ub,threshold,K);
	cout << "finish init" << endl;
	finished = false;
	while(!finished){
		clock_t start = clock();
		finished = query_handle.run_query(result);
		cudaDeviceSynchronize();
		double time = (double)(clock() - start) / CLOCKS_PER_SEC;
		total_time += time;
		cout << "iteration takes time: " << time << endl;
	}

	cout << "finished with total time : " << total_time << endl;

	for(int i = 0; i < result.size(); i++){
		result[i].print_result_entry();
		cout << endl;
	}

	cout << "timing" << endl;
	for(int i = 0; i < exec_time.size(); i++){
		cout << "part " << i << " takes " << exec_time[i] << endl;
	}

	//init_problem("../data/sift_exp_data/","../data/sift_exp_data/query_format1/002_0001.sift",0);
	//cudaProfilerStart();

	//cudaProfilerStop();
	//print_final_result(result1);
	
	//print_final_result(result2);

	
	return 1;
}